Joint density-functional theory for electronic structure of solvated systems 
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We introduce a new form of density functional theory for the ab initio description of electronic 
systems in contact with a molecular liquid environment. This theory rigorously joins an electron 
density-functional for the electrons of a solute with a classical density-functional theory for the liquid 
into a single variational principle for the free energy of the combined system. A simple approximate 
functional predicts, without any fitting of parameters to solvation data, solvation energies as well 
as state-of-the-art quantum-chemical cavity approaches, which require such fitting. 



Ab initio electron density- functional methods have 
proved a computationally efficient and accurate tool for 
the exploration of a wide range of issues in condensed- 
matter physics and chemistry. (See, for instance^i^.) 
However, application of this approach is largely lim- 
ited to the solid-state, gas phase chemistry, or to sur- 
face chemistry in high vacuum environments, with prac- 
tical problems involving liquid environments largely out 
of reach. The reason for this unfortunate limitation is 
that proper description of the effects of water demands 
not only the inclusion of many solvent molecules but 
also thermodynamic sampling of many configurations of 
those molecules so as to properly capture the structure 
and response of the liquid over experimental length and 
time scales. Thus, the applicability of such approaches to 
the vast array of problems involving liquid environments, 
from liquid fuel-cell research to biochemistry, is severely 
curtailed. 

In response, much effort (many thousands of articles 
and a large number of reviews j^i^i^i^ and others, in the 
last decade alone) has been dedicated to the development 
and application of a large number of different "contin- 
uum" solvation models, in which the details of the molec- 
ular aqueous environment are replaced by a continuum 
description. More recently, Dzubiella et ali^ developed a 
continuum model for water where hydrophobic, disper- 
sion and electrostatic energy terms are written as func- 
tionals of the exclusion volume. The plethora of models 
evidences the importance of the problem, but the lack of a 
consensus underscores that no truly satisfactory method 
has been found. 

The weaknesses of the current state of the art in con- 
tinuum solvation approaches arises from their basic struc- 
ture. As Refii and the other reviews cited above describe, 
such approaches generally divide the free energy of sol- 
vation into a number of contributions, typically cavita- 
tion (formation of the solvent-solute interface), disper- 
sion (long-range electron correlations), repulsion (short- 
range electron overlap effects) , and electrostatic (reorien- 
tational and polarization screening in the solvent). Work 
then proceeds to develop models for each of these terms 
separately. The models for these terms generally re- 
quire a cavity shape, which is most usually defined by 
spheres representing atoms or functional groups, where 
the sphere radii are determined by fits of the final re- 
sults to a large empirical database of solvation energies^. 



Finally, to account for nonlinear saturation effects near 
solutes, an intermediate dielectric constant is often used 
in a shell around the solute as a buffer between the solute 
and bulk medium, as in Refi^, for example. Despite the 
successes of this method, the ad hoc separation of physi- 
cal effects (all originating ultimately from the underlying 
quantum and statistical mechanics) and the empirical fit- 
ting to precisely the class of quantities of interest limits 
the reliability of the predictions of these models. This is 
especially true for applications to new classes of chemi- 
cal species or to situations outside of the fitting database, 
such as to study liquid phase surface catalysis. 

The new approach which we propose to pursue here 
suffers none of the aforementioned difficulties. We first 
prove a theorem which shows that the thermodynamics 
of a system and its electrons (solute) in equilibrium with 
a liquid environment (solvent) can be described rigor- 
ously in terms of a joint density- functional theory ( JDFT) 
between the electrons in the system and the molecules 
comprising the solvent. The physics of the equilibrium 
between a solute and a solvent (cavity formation, dielec- 
tric screening, dispersion and repulsion) are then all de- 
termined in a single variational principle. Maintaining 
the first principles nature of density-functional theory, 
this new approach thus requires no artificial separation 
of contributions, no ad hoc definitions of cavity shapes, 
and no empirical fitting of parameters to experimental 
solvation energies. As shown below, even a preliminary 
implementation of this new approach gives results which 
are competitive with state-of-the-art continuum solvation 
models, even without direct fitting to any solvation data. 
This suggests that further refinements will result in a 
new, efficient and predictive approach to electronic struc- 
ture in the presence of liquid environments. 

Joint density-functional theory — A straightforward 
combination of Mermin's non-zero temperature formula- 
tion of density-functional theoryi^ with Capitani et al.^s 
extensions of the zero-temperature theory to include nu- 
clear degrees of freedomii leads to the following, exact 
variational principle for the total thermodynamic free en- 
ergy of an electron-nuclear system in a fixed external elec- 
trostatic potential V{r) 
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min 
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where nt{r) is the thermahy and quantum mechanicaUy 
averaged total single-particle number density of elec- 
trons, Ni{r) is the likewise averaged density of the nu- 
clear species i (of atomic number Z,), and F is a univer- 
sal functional. (Here, as throughout this work, we em- 
ploy atomic units, in which Planck's constant and the 
mass and charge of the electron all have value unity, 
h = nif. = e = 1.) The universality properties of the 
functional F may be seen directly from its construction 
within Levy's constrained search procedurei^, 

F[nt{r),{N.{r)}]^ (2) 
min Tr ( pH -\- ksT phi p] , 

where fcsT is the thermal energy, H represents all in- 
teractions among and the kinetic energy of the electrons 
and nuclei, p is the full quantum-mechanical density ma- 
trix for the electron and nuclear degrees of freedom, and 
the minimization is carried out over only those p which 
lead to the given densities nt{r) and {iVi(r)}. From this 
construction it is clear that F, like H from which it de- 
rives, is independent of the external potential V{r) and 
depends only upon the identities of the nuclear species i 
(and, implicitly, upon the temperature T), as Capitani 
et alfii found previously for the case of T = 0. 

To employ the variational principle Eq. ([T]) in the study 
of a system to be treated explicitly in contact with a sol- 
vent environment, we take the nuclear species i to be 
those comprising the environment (solvent) and the po- 
tential V{r) in Eq. ([T]) to be that arising from the nuclei 
of the explicit system (solute), 



Rj 



(3) 



which we take to sit at fixed locations Rj with atomic 
numbers Zi. Note that, although we employ a Born- 
Oppcnheimer approximation for the nuclei of the explicit 
system, at this stage the treatment of the nuclear species 
of the environment in Eq. ([2]) is fully quantum mechan- 
ical. Thus, Eqs. (|ll2p account for all zero-point motion 
effects associated with lighter nuclear species in the sol- 
vent, such as may be associated with the protons in liq- 
uid water or with the helium atoms in superfluid helium 
when used as a solvent. 

Although Eqs. pi2p give a rigorous continuum treat- 
ment of the environment nuclei, the variational principle 
Eq. ([T]) is ultimately impracticable because it requires ex- 
plicit treatment of all of the electrons, including those as- 
sociated with the environment. We thus "integrate out" 
the electrons associated with the environment by writing 
n((r) = n{r) + ne(r), where ne(r) is the electron density 
associated with the environment and n{r) is the electron 
density associated with the system in contact with that 



environment. We then perform the minimization over 
all allowable ne{r), and finally perform the minimiza- 
tion over all allowable n{r). For this purpose, we define 
the sets of allowable nt{r), ne{r) and n{r) to be all N- 
representable functions satisfying the criteria of Gilbert — 
and integrating to the appropriate number of electrons 
for the respective system. Because all thus defined nt{r) 
can be constructed as the sum of some allowable ne(r) 
and some allowable n(r) and because all such allowable 
ne{r) and n(r) sum to an allowable nt(r), this procedure 
is guaranteed to recover the final free energy in Eq. ([1]). 
Thus, we have 



A - niin \G[n{r),{N,ir)},Vir)] (4) 

(fr V{r)n{r) 

where V{r) is defined above in Eq. ([3]) and where 

G[n{r),{N,{r)},V(r)]= (5) 

min \ F[n{r) + ^^(r), {N,{r)}] 

+ j <fr V{r) Z,N,{r) - n,{r)^ | 

is universal in the sense that its functional form, like F 
from which it derives, depends only on the nature of the 
solvent and, implicitly, the temperature, and that its de- 
pendence on the solute is only through the electrostatic 
potential of the nuclei in V{r). The choice to separate 
the interaction — / <PrV{r)n{r) in Eq. ([4]) from the def- 
inition of G limits the interactions which the unknown 
functional G must describe, easing the task of finding 
good approximations. Note, for instance, that with this 
choice V{r) in Eq. ^ now couples to a neutral charge 
distribution, thereby limiting to the greatest extent pos- 
sible the dependence of G on V{r). 

Eq. ^ gives the exact free energy and exact configura- 
tion of the solvent {iVi(r)}. However, care must be taken 
in the interpretation of the n{r) which yield the mini- 
mum value. The indistinguishability of electrons implies 
that there can be no fundamentally meaningful assign- 
ment of electrons as belonging either to the environment 
or to the system, and thus no exact formulation can give 
a unique result for n{r) without some additional pre- 
scription. Indeed, for the exact {Ni{r)} and any n{r) 
which integrates to the correct number of electrons and 
is everywhere less than the exact solution nt{r) so that 
nf,{r) = nt{r) — n{r) is allowable in the above sense, the 
minimization in Eq. ([S]) will find ne(r) = nt{r)-'n{r) and 
thus ultimately produce the exact value for A. There 
is thus a large set of n(r) which yield the minimum 
value in Eq. and the variational principle embod- 
ied in Eqs. (|4l5p satisfies the fundamental condition of 
not enforcing any particular, arbitrary decomposition of 
the total electron density into solvent and environment 
contributions. 
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In practice, however, we expect approximations to 
Eq. ^ to break the above degeneracy and to pick out a 
unique solution. The standard pseudopotential method, 
which replaces the effects of the nuclei and (relatively) 
inert core electrons of a solid or molecule with an effec- 
tive or "pseudo-" potential^, is in fact an example of an 
approximation which tracks only a portion of the total 
electron density and provides results approaching chemi- 
cal accuracy while suffering no pathologies related to the 
underlying degeneracy of an apportionment of electrons 
between two subsystems. 

The existence and reliability of so called "molecular 
pseudopotential" Hamiltonians^ii^ii^ implies the exis- 
tence of reliable approximations to Eq. ([5]) which pick out 
a unique solution for n{r). Such pseudopotential Hamil- 
tonians replace the effects of the nuclei and electrons of a 
collection of molecules on the electrons of an external sys- 
tem (solute) with an effective potential V^{^.-j(r), which 
depends explicitly on the locations of the molecular nu- 
clei {Ri}. Such Hamiltonians have proved to be quite 
accurate. Using them, for instance, Vaidehi et al. find 
the solvation energy of Li'^ to within 0.6 kcal/mole, and 
Kim, Park and co-workers find results for total energies 
with an accuracy acceptable for the study of the problem 
of an excess electron solvated in water. 

Formulating the exact thermodynamics of such Hamil- 
tonians with the same approach that leads to Eq. ^ 
gives directly the principle in Eq. (|4]), but now with 



G[n{r),{N,{r)},V{T)] 



(6) 



min Tr 

p-,[n(r),{Ndr)}] 



pH{R.^},{Rj} + ksTplnp 



where n{r) represents the electron density associated 
with the solute alone and H^R.j^^jijy represents the in- 
ternal electron gas Hamiltonian for the solvent electrons 
alone, the interaction of these electrons with the molecu- 
lar pseudopotential V{/{.}(r) and a model potential func- 
tion U{{Ri}, {Ri}) describing the interactions among the 
environment molecules and the interaction between those 
molecules and the nuclei of the solute through the electro- 
static potential V{r) defined in Eq. ([3]). Because the elec- 
trons have been already apportioned between the solute 
and the solvent during the construction of {^^j, 
the functional G in Eq. ^ represents an example of an 
approximation to Eq. ^ which is both reliable and free 
of any potentially pathological issues associated with de- 
generate solutions for n{r). 

With the functional dependence of G established in 
Eq. we next separate out known components and 
leave an unknown part to be approximated, 

G[n{r),{N,{r)},V{r)] = AKsHr)] 

+ A,[{N,ir)}] + U[n{r), {7V,(r)}, V{r)], (7) 

where Ai<-5[ri(r)] is standard universal Kohn-Sham 
electron-density functional of the explicit system when 
in isolation, Aiq[{Ni{r)}] is the "classical" density- 
functional for the liquid solvent environment when in iso- 
lation, and U[n{r), {Ni{r)},V{r)], defined formally and 



exactly as the difference between the exact functional and 
the sum of the two former functionals, is a new functional 
describing the coupling between the systems. The new 
functional U[n(r), N(r),V{r)] has the same universality 
properties as the functional G from which it derives. 

Construction of approximate functionals — Working 
with ([7]) requires an approximate functional Aiq[{Ni{r)}] 
for the bulk liquid. For water this is an active area of 
researcbiiiiSiiiS. To describe water in our preliminary 
implementation, we built on the ideas and perspective 
of Sun^^. Wc first imagine minimizing over the proton 
density so that a single field remains, the density N{r) of 
the oxygen nuclei, which one may view as the "molecular 
density" as determined by taking the oxygen nucleus to 
define the location of each molecule. Our version of the 
resulting functional then takes the form 

MN{r)] = (8) 
Ad[N{r)] + [ N{r) [ehs{N{r)) - aiV(r)] d^r 



b / N{r) 



g„{r-r')N{r')(fr' -N{r) 



The first term in Eq. ©, A,d[N{r)] = 
ksT J N{r) {\n{N{r)X^) - l) (fr is the analytically 
exact functional for the ideal gas, where fc^T is the 
thermal energy and A is the thermal de Broglie wave 
length of the solute molecules. In the second term, 
Chs = ksT ((3/2) ((1 - ry)-2 - l) - In (1 - r/)) is the 
Percus-Yevick approximation^^ for the free energy per 
particle of a system of hard spheres of diameter d over 
and above that of the ideal gas, where = {ird^ /6)N{r) 
with d being the sphere diameter (fit to experimental 
parameters below). (In retrospect, this work should have 
employed the Carnahan-Starling approximation, which 
more accurately describes hard spheres and represents 
little or no computational overhead. For the range 
of parameters relevant here, the two approximations 
agree to within about 4%, and so we expect little 
significant change to the final results.) The constant 
a in this second term describes the cohesive tendency 
between molecules that holds the fiuid together. The 
third term is written so that it is non-zero only if 
N(r) is not constant, where 6 is a coupling constant 
and .gcr(f) = exp(— r^/2cr^)/(27rCT^)^/^ is a normahzed 
Gaussian describing the range of non-local behavior. 

The first two terms in Eq. ([8]) capture the properties of 
the bulk fluid. The two parameters a and d in these terms 
were adjusted to reproduce the equilibrium density and 
bulk modulus of water with results summarized in Ta- 
blelH For the final term which describes inhomogeneities, 
the non-local coupling strength b and the range a were 
adjusted to reproduce the macroscopic surface tension of 
water 7 and the approximate location Rh of the "bend" 
(as measured by the point of maximum downward cur- 
vature) of the surface tension versus sphere radius pre- 
diction of the molecular dynamics data of ten Wolde and 
co-workers- Table U summarizes these results as well. 
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parameters 


a d 

(^^) (nm) 
mole"' 


b G 

(in?) 


values 


0.3944 0.2918 


0.1561 0.4388 


properties 


Nb B 
(kg/m3) (GPa) 


7 Rb 
(mJ/m^) (nm) 


this work 
experiment (*) 


998.3 2.184 
998.2 2.190 


71.93 3.70 
72.75 ~4 



TABLE I: Fit parameters from ^ with comparisons between 
model and experimental results (at standard conditions of 
20° C and atmospheric pressure. (*) The value for Rt is 
theoretical^. (See text.) 

We emphasize that although some of the parameter fits 
in this preliminary formulation are empirical, they do not 
involve any direct fit to molecular solvation free energies. 

For the coupling functional U[N{r),n{r),V{r)] in ([7]), 
we proceed by dividing it into two parts: long-range 
dielectric screening /S.U sc[n{r) ,V [r)] capturing the ten- 
dency of the molecules in the liquid to be found in orien- 
tations and polarization states that tend to screen long- 
range electric fields, and a short-range electron-overlap 
interaction MJei[N{r)^n(r)]. The long-range screening 
depends only on the electrostatics of the solute and so 
depends only on its electron density n{r) and the nu- 
clear electrostatic potential V{r). The electron-overlap 
contribution depends upon contact between the solvent 
molecules and the solute electron density and so to some 
approximation depends only upon these two densities. 

The lowest-order form for the electronic coupling be- 
tween the liquid and the solute which is compatible with 
translational and rotational symmetry is then 

AC/e;[7V(r),n(r)] = j (f' R j d\ n{r)Vps{r ~ R)N{R). 

Such a lowest-order coupling is a reasonable starting 
point as the overlap is small. With this form, the convo- 
lution kernel Vps plays the role of the average potential 
which an electron at point r feels from a water molecule 
at point R, similar to the "molecular" pseudopotential of 
the type introduced by Cho and coworkers in the different 
context described abovei^. The main difference between 
this potential and that of Ref.^^ is that in the present, 
preliminary formulation, the pseudopotential contains no 
information about the orientation of the molecules and 
so represents some sort of orientational average. To opti- 
mize numerical convergence, we choose to fit Vps (r) to the 
sum of two origin-centered Gaussians of adjustable width 
and amplitude for a total of four adjustable parameters, 

Vps{r)=Aie~^+A2e~^ 

We then adjusted these parameters to reproduce the 
orientationally averaged interaction of a water molecule 
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FIG. 1: Comparison of energy of interaction of water molecule 
with an atom of neon as a function of neon-oxygen distance: 
results from orientation independent pseudopotential (solid 
curve), orientally averaged ab initio data (centers of error 
bars), typical (la) range of values as a function of orienta- 
tion (range of error bars). 

with an atom of neon as a function of distance. Figure [1] 
summarizes the results. The parameters that we found 
are 

Ai = 0.0765; cri = 2.045; A2 = -0.065; era = 2.165. 

With this simple form, we were able to reproduce 
the average interaction to within 1 millihartree (~ 
0.63 kcal/mole) for all distances beyond 2A. (Smaller dis- 
tances are not very relevant as the interaction becomes 
very repulsive. For comparison, at room temperature, 
ksT = 0.93 mH.) Here, there is no fitting to empirical 
data, and the pseudopotential is truly ab initio in the 
same sense as an electron pseudopotential. As such, Vps 
does not afford an opportunity to fold solvation data into 
the construction of our functional. 

Next, the screening term depends upon the long-range 
electrostatic potential of the solute. We take it to be 
the change in electrostatic energy of this potential in the 
presence of a linear, non-local dielectric function of the 
form 

e(r, r') ^ Sir - r') + '^N{r)f{r - r')N{r'), 

where 5{r) is the Dirac-delta function, /(r) is a short- 
ranged function integrating to one, and Nj, and Xb are 
the density and static dielectric polarizabilities, respec- 
tively, of the bulk liquid. This choice ensures a smooth 
transition from the dielectric constant of the bulk to that 
of vacuum over the length-scale of /(r). The choice to 
include both N{r) and N{r') is motivated both by the 
need for e(r, r') to be symmetric and by the notion that 
the response at point r to a field at point r' depends on 
the density of molecules at both locations. The connect- 
ing function /(r) was chosen to be a Gaussian of width 
(7=2.25 bohr=~1.190 A, somewhat larger than the 0-H 
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FIG. 2: Short-ranged repulsive potential added to prevent 
collapse of liquid density N{r) into the strong electric fields 
within the atomic cores. Once prevention of this artificial 
collapse is achieved, the final results are insensitive to the 
choice of this potential. 



bond distance in water. The motivation for using this 
value of a is that, at lower values, dielectric screening 
at short length-scales is so effective that the system may 
lower its energy by bringing fluid density N{r) into the 
atomic cores and the system thus becomes numerically 
unstable. We stress, however, that once stability was 
achieved, the final results were not sensitive to the choice 
of this parameter (typically 10% variation in the solva- 
tion energy over the range ct = 1.5 bohr to a — 2.5 bohr) 
and that this parameter was in no way adjusted to re- 
produce experimental solvation energies. In the future, 
a direct description of the electric polarization in terms 
of the orientational state of the solvent would capture 
dielectric effects directly and remove the need for con- 
structing such a simple, ad hoc model. 

Finally, in conjunction with a nonlocal dielectric re- 
sponse, a term must be added to the energy functional to 
help prevent the aforementioned numerical issues associ- 
ated with penetration of the solvent density into the cores 
of the atoms. To avoid this, we added a short-ranged re- 
pulsive potential of the form J Vrep{r)N{r) Sr inside the 
atomic cores to prevent the overlap of the solvent with 
the nuclei, where Vrepir) is taken to be a rapidly de- 
caying exponential function (rounded within ^0.25 A of 
the origin) of constant 8.44 A~^, leading to an apparent 
hard wall at thermal energy scales near 1.5 A(Figure[2]). 
Again, once the system is numerically stable, the results 
are not sensitive to the form of the repulsion (typically 
1% changes in solvation energy) so long as the repulsion 
effectively cuts off before the natural point of closest ap- 
proach of the solvent at ~2 A. (See Figure [3] and the 
corresponding discussion below.) 

Results and Conclusions — The meaningful physical 
quantities predicted by joint-density functional theory 
include not only the free energy A but also the liquid 
density N{r). Under certain circumstances, the later is 



accessible directly in experiments. For instance, when 
studying the solvation of a water molecule with liquid 
water, the density N{r) in our formulation gives the den- 
sity of oxygen nuclei given a fixed location for one water 
molecule. As observed by Percus^^, the spherical av- 
erage of N{r) thus corresponds to the oxygen-oxygen 
pair distribution function goo{f) measured in experi- 
ments. Figure [3] shows our results for this quantity. In 
good agreement with both x-ray and neutron scattering 
experiments^^, we find N{r) to be essentially zero un- 
til a radius of ~2.0 A, at which point the density rises 
rapidly, overshooting the bulk density before finally ap- 
proaching it. The experiments do show a much more 
pronounced peak and much more structure in the form 
of oscillations which occur beyond 2.5 A. We believe that 
these discrepancies are due to the over simplified model 
we are using for Aiq[N{r)] in ([8]) and that a better liquid 
density- functional would improve this. 

Generally, the coordination number for a liquid is de- 
fined as the integral of goo from r = to the location 
of the first minimum after the coordination peak. In our 
case there are no such oscillations. However, if we de- 
fine the coordination number as the integral of goo up 
to 3.6 A, the last point before the bulk value of unity is 
obtained, we find a coordination of 5.3, relatively close to 
the value of approximately four measured in liquid water 
and far from the close packed value near twelve typically 
found in simple fluids. Thus, our model, as simple as 
it is, captures enough of the physics of water to reflect 
the hydrogen bonding network which leads to tetrahedral 
coordination. The model also appears to reflect the cor- 
rect energetics and to give a correct ah initio prediction 
of the boundary position of the cavity outside of which 
the fluid is excluded. To our knowledge, this is the first 
accurate determination of such a boundary directly from 
first principles. As solvation energies are known to be 
quite sensitive to the construction of the boundary, this 
success gives a strong advantage to the current approach. 

Finally, Figure [3] summarizes the comparison of the 
free energies predicted by the implementation described 
above with both experiment and the predictions of 
state-of-the-art continuum solvation models. Our joint 
density-functional theory results are clearly superior to 
dielectric-cavity-only calculations and are arguably bet- 
ter than state-of-the-art continuum methods that include 
cavity corrections. It is satisfying to see that, without 
any fitting or ad hoc adjustments, the hydrophobicity of 
methane is predicted correctly. It is also quite encourag- 
ing to see that the value we predict for the most basic 
test, water in water, falls much closer than do the quan- 
tum chemical methods to the correct value for inserting 
a water molecule at a fixed location into liquid water. 
While improvements and further studies are needed, the 
quality of these results without adjustment of parameters 
is very encouraging. 

Although the results of the above preliminary imple- 
mentation have the advantage of involving no direct fit- 
ting to solvation energies and are of comparable quality 
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FIG. 3: Spherical average of N{r) about the oxygen nucleus 
of an explicit water molecule in solution scaled to the bulk liq- 
uid density, corresponding to the correlation function goo{r) 
measured in experiments on liquid water. 



to the state-of-the-art, there are a number of weaknesses 
in the preliminary implementation, namely (1) the lack 
of orientational dependence in the pseudopotential, (2) 
an artificial separation of the dielectric response from 
the internal orientational correlations of the liquid, (3) 
an over simplified, linear description of the dielectric re- 
sponse with an ad hoc nonlocal length scale, (4) the re- 
sulting need to artificially prevent the solvent from over- 
lapping the nuclei. All of these can be addressed by a 
formulation including some information describing orien- 
tational ordering in the fluid, which would allow for the 
use of orientation-dependent pseudopotcntials and could 
be used to more naturally describe the dielectric response 
with its non-local and non-linear effects. Development of 
such an approach is the subject of future work. 
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FIG. 4: Computed (vertical axis) versus experimental (hor- 
izontal axis) solvation energies for water, ethanol, methanol 
and methane (from left to right): perfect agreement (diago- 
nal line), published quantum chemistry values (— for all but 
water, for water) with dielectric contribution only (open 
squares) and including cavitation terms (closed squares) , pre- 
liminary results from joint density functional theory (open 
circles) . 
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